##############################################################
##############################################################
### Replication Materials for
### Stefan Müller and Liam Kneafsey:
### Evidence for the Irrelevance of Irrelevant Events
### Political Science Research and Methods
###
### Please get in touch with the authors if you have any questions: 
### stefan.mueller@ucd.ie
### Note: 0000_readme.pdf contains 
### instructions and details about each script and dataset
##############################################################
##############################################################

## load packages
library(dplyr)      # CRAN v1.0.5
library(ggplot2)    # CRAN v3.3.3
library(estimatr)   # CRAN v0.30.2
library(broom)      # CRAN v0.7.4
library(here)       # CRAN v1.0.1
library(parameters) # CRAN v0.13.0
library(stringr)    # CRAN v1.4.0

here::here()
## alternatively, set working directory manually
## setwd("...")

## load custom ggplot2 theme
source("function_theme_base.R")


## load local election results
dat_local_raw <- readRDS("data_elections_local_full.rds") 

## remove candidates who did not compete in two subsequent elections
dat_local <- filter(dat_local_raw, !is.na(diff_votes))

dat_general_raw <- readRDS("data_elections_general_full.rds")

## remove candidates who did not compete in two subsequent elections
dat_general <- filter(dat_general_raw, !is.na(diff_votes))


## keep only treated incumbents in general elections
dat_general_incumbents_treated <- dat_general %>% 
    filter(winner_loser_before != "Untreated") 

## keep only treated candidates in local elections
dat_local_incumbents_treated <- dat_local %>% 
    filter(winner_loser_before != "Untreated") 


## keep only treated candidates 
dat_general_partyincumbent_treated <- dat_general %>% 
    filter(winner_loser_before != "Untreated") 



## run models for equivalence test and retrieve the coefficient of winner_loser_before

## general elections: effect of win on support for incumbents
lm_cand_general_share <- estimatr::lm_robust(diff_share ~ 
                                                 winner_loser_before * incumbent +
                                                 party_recoded + election_id + 
                                                 county_gaa,
                                             clusters = county_gaa,
                                             se_type = "stata",
                                             data = dat_general_incumbents_treated)

lm_cand_general_share_tidy <- broom::tidy(lm_cand_general_share) %>% 
    filter(term == "winner_loser_beforeWin:incumbentIncumbent (Candidate elected in t-1)")




## general elections: effect of win on support for incumbents
lm_eq_general <- lm(diff_share ~ 
                        winner_loser_before * incumbent +
                        party_recoded + election_id + 
                        county_gaa,
                    data = dat_general_incumbents_treated)


lm_eq_local <- lm(diff_share ~ 
                      winner_loser_before * incumbent + 
                      party_recoded + election_id + 
                      county_gaa,
                  data = dat_local_incumbents_treated)

## relevel variable to ensure that the interaction coefficient reports candidate from government party * winner
dat_general_partyincumbent_treated$incumbency_status_party <- relevel(factor(dat_general_partyincumbent_treated$incumbency_status_party),
                                                                      ref = "Party: Opposition")

lm_eq_partyincumbent <- lm(diff_share ~ 
                               winner_loser_before * incumbency_status_party +
                               party_recoded + election_id + 
                               county_gaa,
                           data = dat_general_partyincumbent_treated)


sd_general <- 0.3 * sd(dat_general_incumbents_treated$diff_share,
                       na.rm = TRUE) 
sd_general


eq_general <- equivalence_test(lm_eq_general,
                               rule = "classic",
                               ci = 0.95,
                               range = c(-sd_general, sd_general)
) %>%
    as.data.frame() %>%
    filter(str_detect(Parameter, "winner_loser_beforeWin:incumbent")) %>% 
    mutate(model = "General Elections: Incumbents x Win")


sd_local <- 0.3 * sd(dat_local_incumbents_treated$diff_share,
                     na.rm = TRUE) 
sd_local

eq_local <- equivalence_test(lm_eq_local,
                             rule = "classic",
                             ci = 0.95,
                             range = c(-sd_local, sd_local),
) %>%
    as.data.frame() %>%
    filter(str_detect(Parameter, "winner_loser_beforeWin:incumbent")) %>% 
    mutate(model = "Local Elections: Incumbents x Win")


sd_partyinc <- 0.3 * sd(dat_general_partyincumbent_treated$diff_share,
                        na.rm = TRUE)
sd_partyinc

eq_partyinc <- equivalence_test(lm_eq_partyincumbent,
                                rule = "classic",
                                ci = 0.95,
                                range = c(-sd_partyinc, sd_partyinc)
) %>%
    as.data.frame() %>%
    filter(str_detect(Parameter, "winner_loser_beforeWin:incumbency_status_partyParty: Government")) %>% 
    mutate(model = "General Elections: Candidates from Government Parties x Win")


eq_combined <- bind_rows(eq_general, 
                         eq_local,
                         eq_partyinc) %>% 
    mutate(estimate = (CI_high + CI_low) / 2)


## Figure A18 ----
ggplot(filter(eq_combined, str_detect(model, "General Elections: Incumbents x Win")),
       aes(x = "", y = estimate, 
           ymin = CI_low, ymax = CI_high)) +
    geom_hline(aes(yintercept = ROPE_low),
               linetype = "dashed", colour = "grey50") +
    geom_hline(aes(yintercept = ROPE_high),
               linetype = "dashed", colour = "grey50") +
    geom_hline(yintercept = 0, colour = "red") +
    geom_point(size = 4) +
    geom_linerange(size = 1.05) +
    facet_wrap(~model) +
    annotate("segment", x = 0.5, xend = 0.5, y = -0.9, yend = -1.55,
             colour = "grey60", size = 0.8, 
             arrow = arrow(angle = 30,
                           length = unit(0.15, "inches"))) +
    annotate("segment", x = 0.5, xend = 0.5, y = 0.9, yend = 1.55,
             colour = "grey60", size = 0.8, 
             arrow = arrow(angle = 30,
                           length = unit(0.15, "inches"))) +
    annotate("text", x = 0.7, y = sd_local - 0.05, hjust = 1,
             colour = "grey50",
             label = "Upper equivalence bound\n(+0.3 SD of DV)") +
    annotate("text", x = 0.7, y = -sd_local  + 0.05, hjust = 0,
             colour = "grey50",
             label = "Lower equivalence bound\n(-0.3 SD of DV)") +
    coord_flip() +
    scale_y_continuous(limits = c(-3, 3), breaks = c(seq(-3, 3, 0.5))) +
    labs(x = NULL, y = "Effect size") +
    theme(axis.ticks.y = element_blank())
ggsave("fig_a18.pdf",
       width = 9, height = 3)


## Figure A19 ----
ggplot(filter(eq_combined, str_detect(model, "Local")),
       aes(x = "", y = estimate, 
           ymin = CI_low, ymax = CI_high)) +
    geom_hline(aes(yintercept = ROPE_low),
               linetype = "dashed", colour = "grey50") +
    geom_hline(aes(yintercept = ROPE_high),
               linetype = "dashed", colour = "grey50") +
    geom_hline(yintercept = 0, colour = "red") +
    geom_point(size = 4) +
    geom_linerange(size = 1.05) +
    facet_wrap(~model) +
    annotate("segment", x = 0.5, xend = 0.5, y = -0.9, yend = -1.4,
             colour = "grey60", size = 0.8, 
             arrow = arrow(angle = 30,
                           length = unit(0.15, "inches"))) +
    annotate("segment", x = 0.5, xend = 0.5, y = 0.9, yend = 1.4,
             colour = "grey60", size = 0.8, 
             arrow = arrow(angle = 30,
                           length = unit(0.15, "inches"))) +
    annotate("text", x = 0.7, y = sd_local - 0.05, hjust = 1,
             colour = "grey50",
             label = "Upper equivalence bound\n(+0.3 SD of DV)") +
    annotate("text", x = 0.7, y = -sd_local  + 0.05, hjust = 0,
             colour = "grey50",
             label = "Lower equivalence bound\n(-0.3 SD of DV)") +
    coord_flip() +
    scale_y_continuous(limits = c(-3, 3), breaks = c(seq(-3, 3, 0.5))) +
    labs(x = NULL, y = "Effect size") +
    theme(axis.ticks.y = element_blank())
ggsave("fig_a19.pdf",
       width = 9, height = 3)


## Figure A20 ----
ggplot(filter(eq_combined, str_detect(model, "Government")),
       aes(x = "", y = estimate, 
           ymin = CI_low, ymax = CI_high)) +
    geom_hline(aes(yintercept = ROPE_low),
               linetype = "dashed", colour = "grey50") +
    geom_hline(aes(yintercept = ROPE_high),
               linetype = "dashed", colour = "grey50") +
    geom_hline(yintercept = 0, colour = "red") +
    geom_point(size = 4) +
    geom_linerange(size = 1.05) +
    facet_wrap(~model) +
    annotate("segment", x = 0.5, xend = 0.5, y = -0.9, yend = -1.55,
             colour = "grey60", size = 0.8, 
             arrow = arrow(angle = 30,
                           length = unit(0.15, "inches"))) +
    annotate("segment", x = 0.5, xend = 0.5, y = 0.9, yend = 1.55,
             colour = "grey60", size = 0.8, 
             arrow = arrow(angle = 30,
                           length = unit(0.15, "inches"))) +
    annotate("text", x = 0.7, y = sd_partyinc - 0.1, hjust = 1,
             colour = "grey50",
             label = "Upper equivalence bound\n(+0.3 SD of DV)") +
    annotate("text", x = 0.7, y = -sd_partyinc + 0.1, hjust = 0,
             colour = "grey50",
             label = "Lower equivalence bound\n(-0.3 SD of DV)") +
    coord_flip() +
    scale_y_continuous(limits = c(-3, 3), breaks = c(seq(-3, 3, 0.5))) +
    labs(x = NULL, y = "Effect size") +
    theme(axis.ticks.y = element_blank())
ggsave("fig_a20.pdf",
       width = 9, height = 3)


## Permutation tests ----


## general elections

#                   # create data frame with one observation per constituency and election 
## and the "treatment"

dat_general_treated <- dat_general %>% 
    mutate(const_election_id = paste(election_id, constituency_name, sep = "_")) %>% 
    filter(winner_loser_before != "Untreated") # remove untreated constituencies

## get a data frame on the level of constituencies that indicates 
## the treatment for this constituency
dat_general_constituencies_treated <- dat_general_treated %>% 
    ungroup() %>% 
    select(election_id, constituency_name, winner_loser_before) %>% 
    unique() %>% 
    arrange(election_id, constituency_name) %>% 
    mutate(const_election_id = paste(election_id, constituency_name, sep = "_"))

nrow(dat_general_constituencies_treated)

## draw 1000 random seeds
set.seed(2356)
seeds <- sample(1:1000, replace = FALSE)

dat_coefficients_permutation_general <- data.frame()

for (i in seeds) {
    
    ## cat("Permutating treatments using seed #", i, "\n")
    set.seed(i)
    
    ## for each election (group_by(election_id)) reshuffle the treatments on the level of
    ## constituencies 
    dat_const_resampled_outcome_general <- dat_general_constituencies_treated %>% 
        group_by(election_id) %>% 
        mutate(winner_loser_before_simulated = winner_loser_before[sample(row_number())], 
               replace = FALSE) %>% 
        ungroup() %>% 
        select(const_election_id, winner_loser_before_simulated) ## select constituency identified and simulated treatment
    
    dat_general_joined <- left_join(dat_general_treated, dat_const_resampled_outcome_general,
                                    by = "const_election_id")
    
    ## run regression model with simulated treatment
    lm_diff_share_sampled_general <- estimatr::lm_robust(
        diff_share ~ 
            winner_loser_before_simulated * incumbent +
            party_recoded + election_id + 
            county_gaa,
        clusters = county_gaa,
        se_type = "stata",
        data = dat_general_joined)
    
    ## tidy regression output
    lm_result_tidy_general <- broom::tidy(lm_diff_share_sampled_general)
    
    ## indicator for seed
    lm_result_tidy_general$seed <- i
    
    ## bind simulations
    dat_coefficients_permutation_general <- bind_rows(lm_result_tidy_general,
                                                      dat_coefficients_permutation_general)
}



## run regression model with real treatment
lm_cand_share_true_general_all <- estimatr::lm_robust(
    diff_share ~ 
        winner_loser_before * incumbent +
        party_recoded + election_id + 
        county_gaa,
    clusters = county_gaa,
    se_type = "stata",
    data = dat_general_treated)

lm_result_tidy_true_general <- broom::tidy(lm_cand_share_true_general_all)


## get coefficient of interaction between simulated treatment and incumbent
dat_coefficients_permutation_general_int <- filter(dat_coefficients_permutation_general,
                                                   term == "winner_loser_before_simulatedWin:incumbentIncumbent (Candidate elected in t-1)") 


## number of simulations (nrow should be 1,000)
nrow(dat_coefficients_permutation_general_int) 

## get coefficient of "true" interaction 
lm_result_tidy_true_general <- filter(lm_result_tidy_true_general, 
                                      term == "winner_loser_beforeWin:incumbentIncumbent (Candidate elected in t-1)")
nrow(lm_result_tidy_true_general)


## Figure A21 ----

#                   # create a plot with the distribution of simulated coefficients
## and the true interaction coefficient
ggplot(dat_coefficients_permutation_general_int, aes(x = estimate)) +
    geom_histogram(colour = "grey50", fill = "grey50") + 
    geom_vline(data = lm_result_tidy_true_general, 
               aes(xintercept = estimate), 
               size = 1.5, 
               colour = "red") +
    geom_vline(xintercept = 0, linetype = "dashed",
               colour = "black") +
    scale_x_continuous(limits = c(-2.3, 2.3)) +
    annotate("text", x = -0.7,
             colour = "white",
             y = 10, label = "Coefficients using\nrandomly permutated treatments") +
    annotate("text", x = lm_result_tidy_true_general$estimate + 0.04,
             hjust = 0,
             colour = "red",
             y = 100, label = "Coefficient using real treatment") +
    labs(x = "Coefficient of interaction term (Win x Candidate elected in t-1)", y = "Count") 
ggsave("fig_a21.pdf",
       width = 9, height = 4)



## repeat the same analysis for local elections

dat_local_treated <- dat_local %>% 
    mutate(const_election_id = paste(election_id, const_link, sep = "_")) %>% 
    filter(winner_loser_before != "Untreated")

dat_local_constituencies_treated <- dat_local_treated %>% 
    ungroup() %>% 
    select(election_id, const_link, winner_loser_before) %>% 
    unique() %>% 
    arrange(election_id, const_link) %>% 
    mutate(const_election_id = paste(election_id, const_link, sep = "_"))

set.seed(2356)
seeds <- sample(1:1000, replace = FALSE)

dat_coefficients_permutation_local <- data.frame()
for (i in seeds) {
    
    ## cat("Permutating treatments using seed #", i, "\n")
    set.seed(i)
    
    ## for each election (group_by(election_id)) reshuffle the treatments on the level of
    ## constituencies 
    dat_const_resampled_outcome_local <- dat_local_constituencies_treated %>% 
        group_by(election_id) %>% 
        mutate(winner_loser_before_simulated = winner_loser_before[sample(row_number())], 
               replace = FALSE) %>% 
        ungroup() %>% 
        select(const_election_id, winner_loser_before_simulated) ## select constituency identified and simulated treatment
    
    dat_local_joined <- left_join(dat_local_treated, dat_const_resampled_outcome_local,
                                  by = "const_election_id")
    
    ## run regression model with simulated treatment
    lm_diff_share_sampled_local <- estimatr::lm_robust(
        diff_share ~ 
            winner_loser_before_simulated * incumbent +
            party_recoded + election_id + 
            county_gaa,
        clusters = county_gaa,
        se_type = "stata", 
        data = dat_local_joined)
    
    ## tidy regression output
    lm_result_tidy_local <- broom::tidy(lm_diff_share_sampled_local)
    
    ## indicator for seed
    lm_result_tidy_local$seed <- i
    
    ## bind rows
    dat_coefficients_permutation_local <- bind_rows(lm_result_tidy_local,
                                                    dat_coefficients_permutation_local)
}


## run regression model with real treatment
lm_cand_share_true_local_all <- estimatr::lm_robust(
    diff_share ~ 
        winner_loser_before * incumbent +
        party_recoded + election_id + 
        county_gaa,
    clusters = county_gaa,
    se_type = "stata",
    data = dat_local_treated)


lm_result_tidy_true_local <- broom::tidy(lm_cand_share_true_local_all)

## get coefficient of simulated interaction term
dat_coefficients_permutation_local_int <- filter(dat_coefficients_permutation_local,
                                                 term == "winner_loser_before_simulatedWin:incumbentIncumbent (Candidate elected in t-1)") 


## number of simulations (nrow should be 1,000)
nrow(dat_coefficients_permutation_local_int)

## get coefficient of "true" interaction 
lm_result_tidy_true_local <- filter(lm_result_tidy_true_local, 
                                    term == "winner_loser_beforeWin:incumbentIncumbent (Candidate elected in t-1)")
nrow(lm_result_tidy_true_local)



## Figure A22 ----

#                   # create a plot with the distribution of simulated coefficients
## and the true interaction coefficient
ggplot(dat_coefficients_permutation_local_int, aes(x = estimate)) +
    geom_histogram(colour = "grey50", fill = "grey50") + 
    geom_vline(data = lm_result_tidy_true_local, 
               aes(xintercept = estimate), 
               size = 1.5, 
               colour = "red") +
    geom_vline(xintercept = 0, colour = "black",
               linetype = "dashed") +
    annotate("text", x = 0.5,
             size = 3,
             colour = "white",
             y = 30, label = "Coefficients\nusing\nrandomly\npermutated\ntreatments") +
    annotate("text", x = lm_result_tidy_true_local$estimate + 0.04,
             hjust = 0,
             colour = "red",
             y = 200, label = "Coefficient using real treatment") +
    scale_x_continuous(limits = c(-2.3, 2.3)) +
    labs(x = "Coefficient of interaction term (Win x Candidate elected in t-1)", y = "Count") 
ggsave("fig_a22.pdf",
       width = 9, height = 4)


## general election, incumbency status of party


## get a data frame on the level of constituencies that indicates 
## the treatment for this constituency

dat_general_partyincumbent_treated <- dat_general_raw %>% 
    filter(winner_loser_before != "Untreated") %>% 
    mutate(const_election_id = paste(election_id, constituency_name, sep = "_"))


dat_general_partyincumbent_treated$incumbency_status_party <- factor(dat_general_partyincumbent_treated$incumbency_status_party,
                                                                     levels = c("Party: Opposition", "Party: Government"))

dat_general_partyincumbent_constituencies_treated <- dat_general_partyincumbent_treated %>% 
    ungroup() %>% 
    select(election_id, constituency_name, winner_loser_before) %>% 
    unique() %>% 
    arrange(election_id, constituency_name) %>% 
    mutate(const_election_id = paste(election_id, constituency_name, sep = "_"))

nrow(dat_general_partyincumbent_constituencies_treated)

## draw 1000 random seeds
set.seed(2356)
seeds <- sample(1:1000, replace = FALSE)

dat_coefficients_permutation_general_partyincumbent <- data.frame()

for (i in seeds) {
    
    ## cat("Permutating treatments using seed #", i, "\n")
    set.seed(i)
    
    ## for each election (group_by(election_id)) reshuffle the treatments on the level of
    ## constituencies 
    dat_const_resampled_outcome_general_partyincumbent <- dat_general_partyincumbent_constituencies_treated %>% 
        group_by(election_id) %>% 
        mutate(winner_loser_before_simulated = winner_loser_before[sample(row_number())], 
               replace = FALSE) %>% 
        ungroup() %>% 
        select(const_election_id, winner_loser_before_simulated) ## select constituency identified and simulated treatment
    
    dat_general_partyincumbent_joined <- left_join(dat_general_partyincumbent_treated, dat_const_resampled_outcome_general_partyincumbent,
                                                   by = "const_election_id")
    
    ## run regression model with simulated treatment
    lm_diff_share_sampled_general_partyincumbent <- estimatr::lm_robust(
        diff_share ~ 
            winner_loser_before_simulated * incumbency_status_party +
            party_recoded + election_id + 
            county_gaa,
        clusters = county_gaa,
        se_type = "stata",
        data = dat_general_partyincumbent_joined)
    
    ## tidy regression output
    lm_result_tidy_general_partyincumbent <- broom::tidy(lm_diff_share_sampled_general_partyincumbent)
    
    ## indicator for seed
    lm_result_tidy_general_partyincumbent$seed <- i
    
    ## bind siimulations
    dat_coefficients_permutation_general_partyincumbent <- bind_rows(lm_result_tidy_general_partyincumbent,
                                                                     dat_coefficients_permutation_general_partyincumbent)
}


## run regression model with real treatment
lm_cand_share_true_general_partyincumbent_all <- estimatr::lm_robust(
    diff_share ~ 
        winner_loser_before * incumbency_status_party +
        party_recoded + election_id + 
        county_gaa,
    clusters = county_gaa,
    se_type = "stata", 
    data = dat_general_partyincumbent_treated)



## get coefficient of interaction between simulated treatment and incumbent
dat_coefficients_permutation_general_partyincumbent_int <- filter(dat_coefficients_permutation_general_partyincumbent,
                                                                  term == "winner_loser_before_simulatedWin:incumbency_status_partyParty: Government") 


## number of simulations (nrow should be 1,000)
nrow(dat_coefficients_permutation_general_partyincumbent_int)

## get coefficient of "true" interaction 
lm_result_tidy_true_general_partyincumbent <- broom::tidy(lm_cand_share_true_general_partyincumbent_all)

lm_result_tidy_true_general_partyincumbent <- filter(lm_result_tidy_true_general_partyincumbent, 
                                                     term == "winner_loser_beforeWin:incumbency_status_partyParty: Government")

## Figure A23 ----

#                   # create a plot with the distribution of simulated coefficients
## and the true interaction coefficient
ggplot(dat_coefficients_permutation_general_partyincumbent_int, aes(x = estimate)) +
    geom_histogram(colour = "grey50", fill = "grey50") + 
    geom_vline(data = lm_result_tidy_true_general_partyincumbent, 
               aes(xintercept = estimate), 
               size = 1.5, 
               colour = "red") +
    geom_vline(xintercept = 0, linetype = "dashed",
               colour = "black") +
    scale_x_continuous(limits = c(-2.3, 2.3)) +
    annotate("text", x = -0.69,
             colour = "white",
             y = 7, label = "Coefficients using\nrandomly permutated treatments") +
    annotate("text", x = lm_result_tidy_true_general_partyincumbent$estimate + 0.04,
             hjust = 0,
             colour = "red",
             y = 100, label = "Coefficient using real treatment") +
    labs(x = "Coefficient of interaction term (Win x Party: Incumbent Government)", y = "Count") 
ggsave("fig_a23.pdf",
       width = 9, height = 4)


## script executed successfully on
date()

sessionInfo()
